1_mw_sleep_and_pubertal_stage
Jackie
06/16/21; 07/06/21
#Load Libraries
library(tidyverse)
library(modelr)
library(foreign)
library(readxl)
library(haven)
library(expss)
library(chron)
library(hms)
library(lubridate)
library(lmerTest)1 Reading in mw and puberty/demo dfs
mw_fp <-
"~/Box/Mooddata_Coordinating/1_Lab_Coordinating/Users/JackieSchwartz/Dissertation/0_MW_Act_Demo_Descriptives/MW_daily_sleep_data_reduced.csv"
# full AHC sample with pubertal and demographic info
pub_demo_act_mw_comb_fp <- "~/Box/Mooddata_Coordinating/1_Lab_Coordinating/Users/JackieSchwartz/Dissertation/0_MW_Act_Demo_Descriptives/dem_and_tanner_mw_act_comb_merge.csv"mw <- read_csv(mw_fp) %>%
mutate(ELS_ID = factor(ELS_ID)) %>%
as_tibble()
pub <- read_csv(pub_demo_act_mw_comb_fp) %>% # remember this dataframe is based on a bunch of scriptss that aligned the tanner and age with the AHC date (all variables - eg, Sex, BMI are based on the AHC timepoint)
mutate(ELS_ID = factor(ELS_ID)) %>%
as_tibble()1.1 filtering for EMA
pub_filt <-
pub %>%
filter(
!is.na(MW_daily_sleep_session_usable)
)1.2 merging dfs
And figuring out if any rows were duplicated by merging
pub_mw <-
left_join(
pub_filt,
mw,
by = "ELS_ID"
) %>%
mutate(
dup = duplicated(ELS_ID)
)
ELSdup <- pub_mw %>% filter(dup == TRUE) # none
write_csv(pub_mw, "pub_mw_merged.csv")1.2.1 functions
getmode <- function(v, na.rm=TRUE) {
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}1.2.2 cronbach’s alpha
library(ltm)
cmep_cron <-
pub_mw %>%
dplyr::select(
starts_with("cmep") &
ends_with("mw"),
-cmep_total_mw
) %>%
drop_na()
cronbach.alpha(cmep_cron) # 0.827##
## Cronbach's alpha for the 'cmep_cron' data-set
##
## Items: 10
## Sample units: 105
## alpha: 0.827# I already computed cronbach's alpha for ASQ in my scoring asq script1.3 Data Viz
1.3.1 distribution of sleep measures
dailysleep_sat_mw <-
pub_mw %>%
ggplot(
aes(x = dailysleep_sat_mean)
) +
geom_histogram(color = "black", alpha = .5, bins =15) +
theme_classic() +
labs(
x = "Averaged Daily Subjective Sleep Satisfaction",
title = "Distribution of Averaged Daily Sleep Satisfaction"
) +
theme(
axis.text = element_text(size = 10, angle = 30, hjust = 1),
axis.title = element_text(size = 10),
plot.title = element_text(size = 10)
)
dailysleep_sat_mw
dailysleep_hrs_mw <-
pub_mw %>%
ggplot(
aes(x = dailysleep_hrs_mean)
) +
geom_histogram(color = "black", alpha = .5, bins =15) +
theme_classic() +
labs(
x = "Averaged Daily Subjective Sleep Duration (hrs)",
title = "Distribution of Averaged Daily Sleep Duration"
) +
theme(
axis.text = element_text(size = 10, angle = 30, hjust = 1),
axis.title = element_text(size = 10),
plot.title = element_text(size = 10)
)
dailysleep_hrs_mw
cmep_mw <-
pub_mw %>%
ggplot(
aes(x = cmep_total_mw)
) +
geom_histogram(color = "black", alpha = .5, bins =15) +
theme_classic() +
labs(
x = "Morningness",
title = "Distribution of Circadian Preference"
) +
theme(
axis.text = element_text(size = 10, angle = 30, hjust = 1),
axis.title = element_text(size = 10),
plot.title = element_text(size = 10)
)
cmep_mw## Warning: Removed 6 rows containing non-finite values (stat_bin).
dailysleep_summary <-
pub_mw %>%
summarise(
dailysleep_sat_average = mean(dailysleep_sat_mean),
dailysleep_sat_disp = sd(dailysleep_sat_mean),
dailysleep_sat_min = min(dailysleep_sat_mean),
dailysleep_sat_max = max(dailysleep_sat_mean),
dailysleep_hrs_average = mean(dailysleep_hrs_mean),
dailysleep_hrs_disp = sd(dailysleep_hrs_mean),
dailysleep_hrs_min = min(dailysleep_hrs_mean),
dailysleep_hrs_max = max(dailysleep_hrs_mean),
cmep_average = mean(cmep_total_mw, na.rm = TRUE),
cmep_disp = sd(cmep_total_mw, na.rm = TRUE),
cmep_min = min(cmep_total_mw, na.rm = TRUE),
cmep_max = max(cmep_total_mw, na.rm = TRUE)
)
dailysleep_summary## # A tibble: 1 x 12
## dailysleep_sat_average dailysleep_sat_disp dailysleep_sat_m… dailysleep_sat_m…
## <dbl> <dbl> <dbl> <dbl>
## 1 60.2 16.8 19.8 100
## # … with 8 more variables: dailysleep_hrs_average <dbl>,
## # dailysleep_hrs_disp <dbl>, dailysleep_hrs_min <dbl>,
## # dailysleep_hrs_max <dbl>, cmep_average <dbl>, cmep_disp <dbl>,
## # cmep_min <dbl>, cmep_max <dbl>getmode(pub_mw$dailysleep_sat_mean) # 78.4## [1] 78.4median(pub_mw$dailysleep_sat_mean) # 62.38462## [1] 62.38462getmode(pub_mw$dailysleep_hrs_mean) # 7## [1] 7median(pub_mw$dailysleep_hrs_mean) # 7.4## [1] 7.4pub_mw_finalcmep <-
pub_mw %>%
drop_na(cmep_total_mw)
getmode(pub_mw_finalcmep$cmep_total_mw) # 23## [1] 23median(pub_mw_finalcmep$cmep_total_mw) # 26## [1] 26library(gridExtra)
library(cowplot)
mwsleep <-
cowplot::plot_grid(
dailysleep_sat_mw,
dailysleep_hrs_mw,
cmep_mw,
labels = c('A', 'B', 'C'),
label_size = 10
)## Warning: Removed 6 rows containing non-finite values (stat_bin).ggsave("averaged_daily_sleep_mw.png", mwsleep, width = 7, height = 5)1.3.2 distribution of most recent tanner score
tanner_mw <-
pub_mw %>%
ggplot(
aes(x = tanner_average_for_mw)
) +
geom_histogram(color = "black", alpha = .5) +
theme_classic() +
labs(
x = "Tanner Stage",
title = "Distribution of Tanner Stage closest to EMA"
) +
theme(
axis.text = element_text(size = 10, angle = 30, hjust = 1),
axis.title = element_text(size = 10),
title = element_text(size = 10),
plot.title = element_text(size = 10)
)
tanner_mw## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
1.3.3 distribution of most recent age
age_mw <-
pub_mw %>%
ggplot(
aes(x = tanner_age_for_mw)
) +
geom_histogram(color = "black", alpha = .5, bins =21) +
theme_classic() +
labs(
x = "Age",
title = "Distribution of Age closest to EMA"
) +
theme(
axis.text = element_text(size = 10, angle = 30, hjust = 1),
axis.title = element_text(size = 10),
title = element_text(size = 10),
plot.title = element_text(size = 10)
) +
scale_x_continuous(
name = "Age",
limits = c(13, 21),
breaks = seq(13, 21, 1)
)
age_mw
1.4 pubertal stage relative to age
within each sex: modeling pubertal stage on age, extracting residuals to get a relative pub stage score
rel_pub_mod <- lm(scale(tanner_average_for_mw) ~ scale(tanner_age_for_mw), data = pub_mw)
summary(rel_pub_mod)##
## Call:
## lm(formula = scale(tanner_average_for_mw) ~ scale(tanner_age_for_mw),
## data = pub_mw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3777 -0.7683 0.1568 0.6621 1.4587
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.899e-16 8.774e-02 0.000 1
## scale(tanner_age_for_mw) 3.915e-01 8.814e-02 4.442 2.15e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9244 on 109 degrees of freedom
## Multiple R-squared: 0.1532, Adjusted R-squared: 0.1455
## F-statistic: 19.73 on 1 and 109 DF, p-value: 2.152e-05pub_mw <-
pub_mw %>%
add_residuals(rel_pub_mod, var = "rel_pub_stg")
relpub_mw <-
pub_mw %>%
ggplot(
aes(x = rel_pub_stg)
) +
geom_histogram(color = "black", alpha = .5) +
theme_classic() +
labs(
x = "Pubertal Stage Relative to Age",
title = "Distribution of Relative Pubertal Stage"
) +
theme(
axis.text = element_text(size = 10, angle = 30, hjust = 1),
axis.title = element_text(size = 10),
title = element_text(size = 10),
plot.title = element_text(size = 10)
)
relpub_mw## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
mwsleep_pub <-
cowplot::plot_grid(
tanner_mw,
age_mw,
relpub_mw,
labels = c('A', 'B', 'C'),
label_size = 10
)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.ggsave("age_and_tanner_mw_dist.png", mwsleep_pub, width = 9, height = 6)options(scipen = 999, digits = 5)1.5 removing outliers in my DV
# outliers: outside 1.5 times the interquartile range above the upper quartile and below the lower quartile (Q1 - 1.5 * IQR or Q3 + 1.5 * IQR).
boxplot.stats(pub_mw$dailysleep_sat_mean)$out # 0 outliers## numeric(0)1.6 centering based on this sample
I standardized tanner stage, age at tanner assessment, and bmi within each sex
pub_mw_cent <-
pub_mw %>%
mutate(
Sex = factor(Sex),
COVID_AHC_MW = factor(COVID_AHC_MW),
mw_dailysleep_sat_z = scale(dailysleep_sat_mean),
mw_dailysleepshrs_z = scale(dailysleep_hrs_mean),
household_income_mw = scale(Household_Income_MW),
sumsev_type_z = scale(sumsev_type),
asq_total_mw_z = scale(asq_total_mw),
rr_z = scale(sleep_resp_rate.x),
) %>%
group_by(Sex) %>%
mutate(
tanner_average_for_mw_z = scale(tanner_average_for_mw), # standardize within each sex
tanner_adrenal_mw_z = scale(tanner_adrenal_mw),
tanner_gonadal_mw_z = scale(tanner_gonadal_mw),
bmi_at_mw_z = scale(BMI_MW),
tanner_age_for_mw_z = scale(tanner_age_for_mw)
) %>%
ungroup()1.7 computing relative pubertal stage
rel_pub_mod <- lm(tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_cent)
summary(rel_pub_mod)##
## Call:
## lm(formula = tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.402 -0.762 0.170 0.656 1.470
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.000000000000000645 0.087659923676139395 0.00 1
## tanner_age_for_mw_z 0.383466939515480276 0.088460487540770100 4.33 0.000033
##
## (Intercept)
## tanner_age_for_mw_z ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.924 on 109 degrees of freedom
## Multiple R-squared: 0.147, Adjusted R-squared: 0.139
## F-statistic: 18.8 on 1 and 109 DF, p-value: 0.0000326pub_mw_cent <-
pub_mw_cent %>%
add_residuals(rel_pub_mod, var = "rel_pub_stg")1.8 associations with sleep satisfaction/restfulness
cor.test(pub_mw_cent$dailysleep_sat_mean, pub_mw_cent$bmi_at_mw_z) # not sig##
## Pearson's product-moment correlation
##
## data: pub_mw_cent$dailysleep_sat_mean and pub_mw_cent$bmi_at_mw_z
## t = -0.243, df = 104, p-value = 0.81
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.21361 0.16769
## sample estimates:
## cor
## -0.023828# r=-0.02439007, p=0.804
cor.test(pub_mw_cent$dailysleep_sat_mean, pub_mw_cent$household_income_mw) # not sig##
## Pearson's product-moment correlation
##
## data: pub_mw_cent$dailysleep_sat_mean and pub_mw_cent$household_income_mw
## t = 1.25, df = 102, p-value = 0.21
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.071223 0.308334
## sample estimates:
## cor
## 0.12305# r=0.1230577, p=0.2133
cor.test(pub_mw_cent$dailysleep_sat_mean, pub_mw_cent$sumsev_type_z) # not sig##
## Pearson's product-moment correlation
##
## data: pub_mw_cent$dailysleep_sat_mean and pub_mw_cent$sumsev_type_z
## t = -0.535, df = 108, p-value = 0.59
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.23635 0.13718
## sample estimates:
## cor
## -0.051385# r=-0.04942569, p=0.6081
cor.test(pub_mw_cent$dailysleep_sat_mean, pub_mw_cent$rr_z) # not sig##
## Pearson's product-moment correlation
##
## data: pub_mw_cent$dailysleep_sat_mean and pub_mw_cent$rr_z
## t = 0.176, df = 109, p-value = 0.86
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.17009 0.20260
## sample estimates:
## cor
## 0.016841# r=0.0123, p=0.9
summary(lm(mw_dailysleep_sat_z ~ Sex, data = pub_mw_cent)) # not sig ##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ Sex, data = pub_mw_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.354 -0.652 0.181 0.676 2.359
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0812 0.1511 0.54 0.59
## Sex2 -0.1345 0.1945 -0.69 0.49
##
## Residual standard error: 1 on 109 degrees of freedom
## Multiple R-squared: 0.00437, Adjusted R-squared: -0.00477
## F-statistic: 0.478 on 1 and 109 DF, p-value: 0.491# B=-0.13231, p=0.498
summary(lm(mw_dailysleep_sat_z ~ COVID_AHC_MW, data = pub_mw_cent)) # not sig##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ COVID_AHC_MW, data = pub_mw_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.460 -0.659 0.186 0.685 2.312
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0527 0.1198 0.44 0.66
## COVID_AHC_MW1 -0.1428 0.1971 -0.72 0.47
##
## Residual standard error: 1 on 109 degrees of freedom
## Multiple R-squared: 0.00479, Adjusted R-squared: -0.00434
## F-statistic: 0.525 on 1 and 109 DF, p-value: 0.47# B=-0.13931, p=0.4811.9 extra libraries
library(gvlma)
library(car)
library(BayesFactor)
library(bayestestR)
library(sjPlot)2 Tanner Average and Sleep Satisfaction Association
2.1 association between relative pubertal stage score and sleep sat
relpub <- lm(mw_dailysleep_sat_z ~ rel_pub_stg, data = pub_mw_cent)
summary(gvlma(relpub))##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ rel_pub_stg, data = pub_mw_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.392 -0.646 0.141 0.738 2.039
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.000000000000000164 0.092254949787475546 0.00 1.0000
## rel_pub_stg -0.274896434084641206 0.100803436898558349 -2.73 0.0074 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.972 on 109 degrees of freedom
## Multiple R-squared: 0.0639, Adjusted R-squared: 0.0553
## F-statistic: 7.44 on 1 and 109 DF, p-value: 0.00745
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = relpub)
##
## Value p-value Decision
## Global Stat 5.628 0.229 Assumptions acceptable.
## Skewness 2.180 0.140 Assumptions acceptable.
## Kurtosis 1.255 0.263 Assumptions acceptable.
## Link Function 2.038 0.153 Assumptions acceptable.
## Heteroscedasticity 0.156 0.693 Assumptions acceptable.summary(relpub)##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ rel_pub_stg, data = pub_mw_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.392 -0.646 0.141 0.738 2.039
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.000000000000000164 0.092254949787475546 0.00 1.0000
## rel_pub_stg -0.274896434084641206 0.100803436898558349 -2.73 0.0074 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.972 on 109 degrees of freedom
## Multiple R-squared: 0.0639, Adjusted R-squared: 0.0553
## F-statistic: 7.44 on 1 and 109 DF, p-value: 0.00745plot_model(relpub, type = 'diag')## [[1]]## `geom_smooth()` using formula 'y ~ x'
##
## [[2]]
##
## [[3]]## `geom_smooth()` using formula 'y ~ x'
# negative association between relative pubertal stage and sleep sat
# (B=-0.2758164704735612816, p=0..0725)
# Overall R2=0.05571, p=0.007246
relpub_int <- summary(relpub)$coefficients[1,1]
relpub_slp <- summary(relpub)$coefficients[2,1]
# CI
ci_relpub = Boot(relpub, coef,
method='case', R=1000)
boot_confint_relpub <- confint(ci_relpub, type='norm')
# 2.5 % 97.5 %
# (Intercept) -0.181 0.1931
# rel_pub_stg -0.478 -0.0724
# Bayes Factor
bayes_relpub <- regressionBF(
formula = mw_dailysleep_sat_z ~ rel_pub_stg, data = pub_mw_cent
) # 5.31## Warning: data coerced from tibble to data frame# plotting
p1 <-
pub_mw_cent %>%
ggplot(
aes(x=rel_pub_stg, y=mw_dailysleep_sat_z)
) +
geom_point(size = 2, alpha = .5) +
geom_abline(
intercept = relpub_int,
slope = relpub_slp,
size=2
) +
labs(
x = "Relative Pubertal Stage",
y = "Averaged Daily Sleep Satisfaction (z-scored)",
title = "Relative Pubertal Stage and Averaged Daily Sleep Satisfaction"
) +
theme_classic() +
theme(plot.title = element_text(size = 9),
axis.text.x = element_text(size = 9),
axis.text.y = element_text(size = 9),
axis.title.x = element_text(size = 9),
axis.title.y = element_text(size = 9))
p1
2.2 separating out tanner and age
# separating out tanner and age
sleep_relpub_mod <- lm(mw_dailysleep_sat_z ~ tanner_average_for_mw_z + tanner_age_for_mw_z, data = pub_mw_cent)
summary(gvlma(sleep_relpub_mod))##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_average_for_mw_z +
## tanner_age_for_mw_z, data = pub_mw_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.357 -0.568 0.154 0.730 1.878
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.0000000000000000183 0.0922621514411763183 0.00
## tanner_average_for_mw_z -0.2748964340846413723 0.1008113058686901792 -2.73
## tanner_age_for_mw_z 0.0131044654376880950 0.1008113058686901375 0.13
## Pr(>|t|)
## (Intercept) 1.0000
## tanner_average_for_mw_z 0.0075 **
## tanner_age_for_mw_z 0.8968
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.972 on 108 degrees of freedom
## Multiple R-squared: 0.0723, Adjusted R-squared: 0.0551
## F-statistic: 4.21 on 2 and 108 DF, p-value: 0.0174
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = sleep_relpub_mod)
##
## Value p-value Decision
## Global Stat 6.13 0.1895 Assumptions acceptable.
## Skewness 2.99 0.0837 Assumptions acceptable.
## Kurtosis 1.34 0.2463 Assumptions acceptable.
## Link Function 1.69 0.1941 Assumptions acceptable.
## Heteroscedasticity 0.11 0.7399 Assumptions acceptable.summary(sleep_relpub_mod) ##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_average_for_mw_z +
## tanner_age_for_mw_z, data = pub_mw_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.357 -0.568 0.154 0.730 1.878
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 0.0000000000000000183 0.0922621514411763183 0.00
## tanner_average_for_mw_z -0.2748964340846413723 0.1008113058686901792 -2.73
## tanner_age_for_mw_z 0.0131044654376880950 0.1008113058686901375 0.13
## Pr(>|t|)
## (Intercept) 1.0000
## tanner_average_for_mw_z 0.0075 **
## tanner_age_for_mw_z 0.8968
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.972 on 108 degrees of freedom
## Multiple R-squared: 0.0723, Adjusted R-squared: 0.0551
## F-statistic: 4.21 on 2 and 108 DF, p-value: 0.0174# still a negative association between tanner stage and sleep satisfaction
# tanner: (B=-0.27581647047356105951, p=0.00726); age: (B=0.01367484884277658450, p=0.89233), overall model Rsq = 0.05553, p=0.01698
sleep_relpub_mod_int <- summary(sleep_relpub_mod)$coefficients[1,1]
sleep_relpub_mod_slp <- summary(sleep_relpub_mod)$coefficients[2,1]
# CI
ci_tanage = Boot(sleep_relpub_mod, coef,
method='case', R=1000)
boot_confint_tanage <- confint(ci_tanage, type='norm')
# 2.5 % 97.5 %
# tanner_average_for_mw_z -0.474 -0.0785
# tanner_age_for_mw_z -0.165 0.2026
# Bayes Factor
bayes_tanage <- regressionBF(
formula = mw_dailysleep_sat_z ~ tanner_average_for_mw_z + tanner_age_for_mw_z, data = pub_mw_cent
) # tanner bf=8.3, age = 0.31## Warning: data coerced from tibble to data frame# plotting
p2 <-
pub_mw_cent %>%
ggplot(
aes(x=tanner_average_for_mw_z, y=mw_dailysleep_sat_z)
) +
geom_point(size = 2, alpha = .5) +
geom_abline(
intercept = sleep_relpub_mod_int,
slope = sleep_relpub_mod_slp,
size=2
) +
labs(
x = "Pubertal Stage (z-scored)",
y = "Averaged Daily Sleep Satisfaction (z-scored)",
title = "Pubertal Stage and Averaged Daily Sleep Satisfaction"
) +
theme_classic() +
theme(plot.title = element_text(size = 9),
axis.text.x = element_text(size = 9),
axis.text.y = element_text(size = 9),
axis.title.x = element_text(size = 9),
axis.title.y = element_text(size = 9))
p2
2.3 gonadal vs. adrenal?
# separating out adrenal and gonadal effects
ad <- lm(mw_dailysleep_sat_z ~ tanner_adrenal_mw_z, data = pub_mw_cent)
summary(gvlma(ad)) # all checks out##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_adrenal_mw_z, data = pub_mw_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.332 -0.638 0.162 0.727 1.804
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.000000000000000317 0.092481486381312636 0.00
## tanner_adrenal_mw_z -0.244562264138701996 0.093326083696018816 -2.62
## Pr(>|t|)
## (Intercept) 1.00
## tanner_adrenal_mw_z 0.01 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.974 on 109 degrees of freedom
## Multiple R-squared: 0.0593, Adjusted R-squared: 0.0506
## F-statistic: 6.87 on 1 and 109 DF, p-value: 0.01
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = ad)
##
## Value p-value Decision
## Global Stat 4.0406 0.401 Assumptions acceptable.
## Skewness 2.5782 0.108 Assumptions acceptable.
## Kurtosis 1.3711 0.242 Assumptions acceptable.
## Link Function 0.0126 0.911 Assumptions acceptable.
## Heteroscedasticity 0.0787 0.779 Assumptions acceptable.summary(ad) ##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_adrenal_mw_z, data = pub_mw_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.332 -0.638 0.162 0.727 1.804
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.000000000000000317 0.092481486381312636 0.00
## tanner_adrenal_mw_z -0.244562264138701996 0.093326083696018816 -2.62
## Pr(>|t|)
## (Intercept) 1.00
## tanner_adrenal_mw_z 0.01 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.974 on 109 degrees of freedom
## Multiple R-squared: 0.0593, Adjusted R-squared: 0.0506
## F-statistic: 6.87 on 1 and 109 DF, p-value: 0.01# B = -0.244562264138701996, p=.01
r_ad <- cor.test(pub_mw_cent$mw_dailysleep_sat_z, pub_mw_cent$tanner_adrenal_mw_z)
r_adcoef <- r_ad$estimate
ci_ad = Boot(ad, coef,
method='case', R=1000)
boot_confint_ad <- confint(ci_ad, type='norm')
# 2.5 % 97.5 %
# tanner_adrenal_mw_z -0.45179 -0.044891
bayes_ad <- regressionBF(
formula = mw_dailysleep_sat_z ~ tanner_adrenal_mw_z, data = pub_mw_cent
) # 4.1549## Warning: data coerced from tibble to data framego <- lm(mw_dailysleep_sat_z ~ tanner_gonadal_mw_z, data = pub_mw_cent)
summary(gvlma(go)) # all checks out##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_gonadal_mw_z, data = pub_mw_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5369 -0.6292 0.0955 0.7405 2.2509
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.000000000000000184 0.093372561247362179 0.00
## tanner_gonadal_mw_z -0.203538707386455686 0.094225296400986397 -2.16
## Pr(>|t|)
## (Intercept) 1.000
## tanner_gonadal_mw_z 0.033 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.984 on 109 degrees of freedom
## Multiple R-squared: 0.0411, Adjusted R-squared: 0.0323
## F-statistic: 4.67 on 1 and 109 DF, p-value: 0.033
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = go)
##
## Value p-value Decision
## Global Stat 4.316 0.3649 Assumptions acceptable.
## Skewness 3.479 0.0621 Assumptions acceptable.
## Kurtosis 0.472 0.4922 Assumptions acceptable.
## Link Function 0.223 0.6368 Assumptions acceptable.
## Heteroscedasticity 0.143 0.7057 Assumptions acceptable.summary(go) ##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_gonadal_mw_z, data = pub_mw_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5369 -0.6292 0.0955 0.7405 2.2509
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.000000000000000184 0.093372561247362179 0.00
## tanner_gonadal_mw_z -0.203538707386455686 0.094225296400986397 -2.16
## Pr(>|t|)
## (Intercept) 1.000
## tanner_gonadal_mw_z 0.033 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.984 on 109 degrees of freedom
## Multiple R-squared: 0.0411, Adjusted R-squared: 0.0323
## F-statistic: 4.67 on 1 and 109 DF, p-value: 0.033# B = -0.203538707386455686, p=0.033
r_go <- cor.test(pub_mw_cent$mw_dailysleep_sat_z, pub_mw_cent$tanner_gonadal_mw_z)
r_gocoef <- r_go$estimate
ci_go = Boot(go, coef,
method='case', R=1000)
boot_confint_go <- confint(ci_go, type='norm')
# 2.5 % 97.5 %
# tanner_gonadal_mw_z -0.39701 -0.012505
bayes_go <- regressionBF(
formula = mw_dailysleep_sat_z ~ tanner_gonadal_mw_z, data = pub_mw_cent
) # 1.5955## Warning: data coerced from tibble to data framep <- c(0.01, 0.033)
p.adjust(p, method = "fdr", n=length(p)) # ad: 0.020 go: 0.033## [1] 0.020 0.033library(cocor)
go_ad <- cor.test(pub_mw_cent$tanner_adrenal_mw_z, pub_mw_cent$tanner_gonadal_mw_z)
goadcoef <- go_ad$estimate
cocor.dep.groups.overlap(r_adcoef, r_gocoef, goadcoef, 111, alternative = "two.sided", test = "all", alpha = 0.05, conf.level = 0.95, null.value = 0, data.name = "pub_mw_cent", var.labels = c("mw_dailysleep_sat_z", "tanner_adrenal_mw_z", "tanner_gonadal_mw_z"), return.htest = FALSE)##
## Results of a comparison of two overlapping correlations based on dependent groups
##
## Comparison between r.jk (mw_dailysleep_sat_z, tanner_adrenal_mw_z) = -0.2434 and r.jh (mw_dailysleep_sat_z, tanner_gonadal_mw_z) = -0.2026
## Difference: r.jk - r.jh = -0.0408
## Related correlation: r.kh = 0.3745
## Data: pub_mw_cent: j = mw_dailysleep_sat_z, k = tanner_adrenal_mw_z, h = tanner_gonadal_mw_z
## Group size: n = 111
## Null hypothesis: r.jk is equal to r.jh
## Alternative hypothesis: r.jk is not equal to r.jh (two-sided)
## Alpha: 0.05
##
## pearson1898: Pearson and Filon's z (1898)
## z = -0.3981, p-value = 0.6906
## Null hypothesis retained
##
## hotelling1940: Hotelling's t (1940)
## t = -0.3942, df = 108, p-value = 0.6942
## Null hypothesis retained
##
## williams1959: Williams' t (1959)
## t = -0.3928, df = 108, p-value = 0.6953
## Null hypothesis retained
##
## olkin1967: Olkin's z (1967)
## z = -0.3981, p-value = 0.6906
## Null hypothesis retained
##
## dunn1969: Dunn and Clark's z (1969)
## z = -0.3926, p-value = 0.6946
## Null hypothesis retained
##
## hendrickson1970: Hendrickson, Stanley, and Hills' (1970) modification of Williams' t (1959)
## t = -0.3942, df = 108, p-value = 0.6942
## Null hypothesis retained
##
## steiger1980: Steiger's (1980) modification of Dunn and Clark's z (1969) using average correlations
## z = -0.3925, p-value = 0.6947
## Null hypothesis retained
##
## meng1992: Meng, Rosenthal, and Rubin's z (1992)
## z = -0.3925, p-value = 0.6947
## Null hypothesis retained
## 95% confidence interval for r.jk - r.jh: -0.2576 0.1717
## Null hypothesis retained (Interval includes 0)
##
## hittner2003: Hittner, May, and Silver's (2003) modification of Dunn and Clark's z (1969) using a backtransformed average Fisher's (1921) Z procedure
## z = -0.3925, p-value = 0.6947
## Null hypothesis retained
##
## zou2007: Zou's (2007) confidence interval
## 95% confidence interval for r.jk - r.jh: -0.2425 0.1619
## Null hypothesis retained (Interval includes 0)ad_int <- summary(ad)$coefficients[1,1]
ad_slp <- summary(ad)$coefficients[2,1]
go_int <- summary(go)$coefficients[1,1]
go_slp <- summary(go)$coefficients[2,1]library(cowplot)
sleepsatplot <-
cowplot::plot_grid(p1,p2, labels = c('A', 'B'), label_size = 8)
sleepsatplot
ggsave("tanner_and_averaged_daily_sleep_sat.png", sleepsatplot, width = 8, height = 5)2.3.0.1 controlling for adolescent stress
cor.test(pub_mw$dailysleep_sat_mean, pub_mw$asq_total_mw) # sig##
## Pearson's product-moment correlation
##
## data: pub_mw$dailysleep_sat_mean and pub_mw$asq_total_mw
## t = -2.69, df = 92, p-value = 0.0086
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.447881 -0.071006
## sample estimates:
## cor
## -0.26974# r=-0.26974, p=0.0086
pub_mw_asq <-
pub_mw %>%
drop_na(asq_total_mw)
pub_mw_asq_cent <-
pub_mw_asq %>%
mutate(
Sex = factor(Sex),
COVID_AHC_MW = factor(COVID_AHC_MW),
mw_dailysleep_sat_z = scale(dailysleep_sat_mean),
mw_dailysleepshrs_z = scale(dailysleep_hrs_mean),
household_income_mw = scale(Household_Income_MW),
sumsev_type_z = scale(sumsev_type),
asq_total_mw_z = scale(asq_total_mw)
) %>%
group_by(Sex) %>%
mutate(
tanner_average_for_mw_z = scale(tanner_average_for_mw), # standardize within each sex
tanner_adrenal_mw_z = scale(tanner_adrenal_mw),
tanner_gonadal_mw_z = scale(tanner_gonadal_mw),
bmi_at_mw_z = scale(BMI_MW),
tanner_age_for_mw_z = scale(tanner_age_for_mw)
) %>%
ungroup()
rel_pub_mod_asq <- lm(tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_asq_cent)
summary(rel_pub_mod_asq)##
## Call:
## lm(formula = tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_asq_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.457 -0.799 0.164 0.711 1.480
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.000000000000000926 0.100201240100561639 0.00
## tanner_age_for_mw_z 0.237092554368757147 0.101284528232989862 2.34
## Pr(>|t|)
## (Intercept) 1.000
## tanner_age_for_mw_z 0.021 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.971 on 92 degrees of freedom
## Multiple R-squared: 0.0562, Adjusted R-squared: 0.046
## F-statistic: 5.48 on 1 and 92 DF, p-value: 0.0214pub_mw_asq_cent <-
pub_mw_asq_cent %>%
add_residuals(rel_pub_mod_asq, var = "rel_pub_stg_asq")
relpubasq <- lm(mw_dailysleep_sat_z ~ rel_pub_stg_asq + asq_total_mw_z, data = pub_mw_asq_cent)
summary(gvlma(relpubasq))##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ rel_pub_stg_asq + asq_total_mw_z,
## data = pub_mw_asq_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4275 -0.6616 0.0155 0.7691 1.8029
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0000000000000000407 0.0975047589756031230 0.00 1.000
## rel_pub_stg_asq -0.2396419814193159747 0.1022501987730255052 -2.34 0.021
## asq_total_mw_z -0.2408587057896816874 0.0987992444877573567 -2.44 0.017
##
## (Intercept)
## rel_pub_stg_asq *
## asq_total_mw_z *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.945 on 91 degrees of freedom
## Multiple R-squared: 0.126, Adjusted R-squared: 0.106
## F-statistic: 6.53 on 2 and 91 DF, p-value: 0.00223
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = relpubasq)
##
## Value p-value Decision
## Global Stat 6.6944 0.1529 Assumptions acceptable.
## Skewness 2.3952 0.1217 Assumptions acceptable.
## Kurtosis 1.0850 0.2976 Assumptions acceptable.
## Link Function 3.1532 0.0758 Assumptions acceptable.
## Heteroscedasticity 0.0609 0.8050 Assumptions acceptable.summary(relpubasq) # relpub: B=-0.23964, p=0.021; asq: B=-0.2408587, p=0.017##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ rel_pub_stg_asq + asq_total_mw_z,
## data = pub_mw_asq_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4275 -0.6616 0.0155 0.7691 1.8029
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0000000000000000407 0.0975047589756031230 0.00 1.000
## rel_pub_stg_asq -0.2396419814193159747 0.1022501987730255052 -2.34 0.021
## asq_total_mw_z -0.2408587057896816874 0.0987992444877573567 -2.44 0.017
##
## (Intercept)
## rel_pub_stg_asq *
## asq_total_mw_z *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.945 on 91 degrees of freedom
## Multiple R-squared: 0.126, Adjusted R-squared: 0.106
## F-statistic: 6.53 on 2 and 91 DF, p-value: 0.00223# R2=0.106, p=0.00223
# CI
ci_relpub = Boot(relpubasq, coef,
method='case', R=1000)
boot_confint_relpub <- confint(ci_relpub, type='norm')
# 2.5 % 97.5 %
# (Intercept) -0.189 0.1917
# rel_pub_stg_asq -0.44263 -0.044854
# asq_total_mw_z -0.42312 -0.069426
# Bayes Factor
bayes_relpub <- regressionBF(
formula = mw_dailysleep_sat_z ~ rel_pub_stg_asq + asq_total_mw_z, data = pub_mw_asq_cent
) # pub: 4.0458 ; asq: 4.9157## Warning: data coerced from tibble to data frame# separating out tanner and age
sleep_relpub_mod <- lm(mw_dailysleep_sat_z ~ tanner_average_for_mw_z + tanner_age_for_mw_z + asq_total_mw_z, data = pub_mw_asq_cent)
summary(gvlma(sleep_relpub_mod))##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_average_for_mw_z +
## tanner_age_for_mw_z + asq_total_mw_z, data = pub_mw_asq_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2963 -0.6261 0.0922 0.7530 1.6662
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.000000000000000394 0.096947253137657730 0.00
## tanner_average_for_mw_z -0.240716832849127732 0.101668331924477440 -2.37
## tanner_age_for_mw_z -0.083468484519627156 0.100910067222319974 -0.83
## asq_total_mw_z -0.232532749472142142 0.098406334375013213 -2.36
## Pr(>|t|)
## (Intercept) 1.00
## tanner_average_for_mw_z 0.02 *
## tanner_age_for_mw_z 0.41
## asq_total_mw_z 0.02 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.94 on 90 degrees of freedom
## Multiple R-squared: 0.145, Adjusted R-squared: 0.117
## F-statistic: 5.09 on 3 and 90 DF, p-value: 0.00267
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = sleep_relpub_mod)
##
## Value p-value Decision
## Global Stat 6.955 0.1383 Assumptions acceptable.
## Skewness 2.543 0.1108 Assumptions acceptable.
## Kurtosis 1.343 0.2465 Assumptions acceptable.
## Link Function 2.968 0.0849 Assumptions acceptable.
## Heteroscedasticity 0.101 0.7505 Assumptions acceptable.summary(sleep_relpub_mod) ##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_average_for_mw_z +
## tanner_age_for_mw_z + asq_total_mw_z, data = pub_mw_asq_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2963 -0.6261 0.0922 0.7530 1.6662
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.000000000000000394 0.096947253137657730 0.00
## tanner_average_for_mw_z -0.240716832849127732 0.101668331924477440 -2.37
## tanner_age_for_mw_z -0.083468484519627156 0.100910067222319974 -0.83
## asq_total_mw_z -0.232532749472142142 0.098406334375013213 -2.36
## Pr(>|t|)
## (Intercept) 1.00
## tanner_average_for_mw_z 0.02 *
## tanner_age_for_mw_z 0.41
## asq_total_mw_z 0.02 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.94 on 90 degrees of freedom
## Multiple R-squared: 0.145, Adjusted R-squared: 0.117
## F-statistic: 5.09 on 3 and 90 DF, p-value: 0.00267# CI
ci_tanage = Boot(sleep_relpub_mod, coef,
method='case', R=1000)
boot_confint_tanage <- confint(ci_tanage, type='norm')
# 2.5 % 97.5 %
# tanner_average_for_mw_z -0.42919 -0.055928
# tanner_age_for_mw_z -0.26290 0.091902
# asq_total_mw_z -0.42506 -0.061145
# Bayes Factor
bayes_tanage <- regressionBF(
formula = mw_dailysleep_sat_z ~ tanner_average_for_mw_z + tanner_age_for_mw_z + asq_total_mw_z, data = pub_mw_asq_cent
) # BF tanner: 8.3485, BF age: 0.5743, BF asq: 4.9157## Warning: data coerced from tibble to data frame# separating out adrenal and gonadal effects
ad <- lm(mw_dailysleep_sat_z ~ tanner_adrenal_mw_z + asq_total_mw_z, data = pub_mw_asq_cent)
summary(gvlma(ad)) # all checks out##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_adrenal_mw_z + asq_total_mw_z,
## data = pub_mw_asq_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2272 -0.6091 0.0557 0.7760 1.7358
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.00000000000000023 0.09810492147505201 0.00 1.000
## tanner_adrenal_mw_z -0.21191364412957212 0.10200592622229875 -2.08 0.041
## asq_total_mw_z -0.22035053767793727 0.10145602503988681 -2.17 0.032
##
## (Intercept)
## tanner_adrenal_mw_z *
## asq_total_mw_z *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.951 on 91 degrees of freedom
## Multiple R-squared: 0.115, Adjusted R-squared: 0.0953
## F-statistic: 5.9 on 2 and 91 DF, p-value: 0.0039
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = ad)
##
## Value p-value Decision
## Global Stat 5.1327 0.274 Assumptions acceptable.
## Skewness 2.2440 0.134 Assumptions acceptable.
## Kurtosis 1.3268 0.249 Assumptions acceptable.
## Link Function 1.4950 0.221 Assumptions acceptable.
## Heteroscedasticity 0.0669 0.796 Assumptions acceptable.summary(ad) ##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_adrenal_mw_z + asq_total_mw_z,
## data = pub_mw_asq_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2272 -0.6091 0.0557 0.7760 1.7358
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.00000000000000023 0.09810492147505201 0.00 1.000
## tanner_adrenal_mw_z -0.21191364412957212 0.10200592622229875 -2.08 0.041
## asq_total_mw_z -0.22035053767793727 0.10145602503988681 -2.17 0.032
##
## (Intercept)
## tanner_adrenal_mw_z *
## asq_total_mw_z *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.951 on 91 degrees of freedom
## Multiple R-squared: 0.115, Adjusted R-squared: 0.0953
## F-statistic: 5.9 on 2 and 91 DF, p-value: 0.0039ci_ad = Boot(ad, coef,
method='case', R=1000)
boot_confint_ad <- confint(ci_ad, type='norm')
# 2.5 % 97.5 %
# tanner_adrenal_mw_z -0.43148 0.0080374
# asq_total_mw_z -0.41442 -0.0260106
bayes_ad <- regressionBF(
formula = mw_dailysleep_sat_z ~ tanner_adrenal_mw_z + asq_total_mw_z, data = pub_mw_asq_cent
) # 4.1235## Warning: data coerced from tibble to data framego <- lm(mw_dailysleep_sat_z ~ tanner_gonadal_mw_z+ asq_total_mw_z, data = pub_mw_asq_cent)
summary(gvlma(go)) # all checks out##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_gonadal_mw_z + asq_total_mw_z,
## data = pub_mw_asq_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.590 -0.666 0.113 0.732 1.929
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.000000000000000114 0.097906561765181491 0.00
## tanner_gonadal_mw_z -0.214632394727764436 0.098976537644006130 -2.17
## asq_total_mw_z -0.272996146822738817 0.098442967516297303 -2.77
## Pr(>|t|)
## (Intercept) 1.0000
## tanner_gonadal_mw_z 0.0327 *
## asq_total_mw_z 0.0067 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.949 on 91 degrees of freedom
## Multiple R-squared: 0.118, Adjusted R-squared: 0.0989
## F-statistic: 6.11 on 2 and 91 DF, p-value: 0.00325
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = go)
##
## Value p-value Decision
## Global Stat 7.6695 0.1045 Assumptions acceptable.
## Skewness 3.3913 0.0655 Assumptions acceptable.
## Kurtosis 0.3716 0.5421 Assumptions acceptable.
## Link Function 3.8279 0.0504 Assumptions acceptable.
## Heteroscedasticity 0.0788 0.7790 Assumptions acceptable.summary(go) ##
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_gonadal_mw_z + asq_total_mw_z,
## data = pub_mw_asq_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.590 -0.666 0.113 0.732 1.929
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.000000000000000114 0.097906561765181491 0.00
## tanner_gonadal_mw_z -0.214632394727764436 0.098976537644006130 -2.17
## asq_total_mw_z -0.272996146822738817 0.098442967516297303 -2.77
## Pr(>|t|)
## (Intercept) 1.0000
## tanner_gonadal_mw_z 0.0327 *
## asq_total_mw_z 0.0067 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.949 on 91 degrees of freedom
## Multiple R-squared: 0.118, Adjusted R-squared: 0.0989
## F-statistic: 6.11 on 2 and 91 DF, p-value: 0.00325ci_go = Boot(go, coef,
method='case', R=1000)
boot_confint_go <- confint(ci_go, type='norm')
# 2.5 % 97.5 %
# tanner_gonadal_mw_z -0.41192 -0.015542
# asq_total_mw_z -0.46247 -0.081226
bayes_go <- regressionBF(
formula = mw_dailysleep_sat_z ~ tanner_gonadal_mw_z+ asq_total_mw_z, data = pub_mw_asq_cent
) # pub: 1.3659, asq: 4.9157## Warning: data coerced from tibble to data frame3 association between sleep satisfaction and sleep duration
cor.test(pub_mw$dailysleep_sat_mean, pub_mw$dailysleep_hrs_mean)##
## Pearson's product-moment correlation
##
## data: pub_mw$dailysleep_sat_mean and pub_mw$dailysleep_hrs_mean
## t = 5.14, df = 109, p-value = 0.0000012
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.27830 0.58038
## sample estimates:
## cor
## 0.44177# r = 0.446, p<0.001
summary(pub_mw$dailysleep_hrs_mean)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.14 6.87 7.40 7.33 8.00 9.00sd(pub_mw$dailysleep_hrs_mean)## [1] 0.85108boxplot.stats(pub_mw$dailysleep_hrs_mean)$out # 1 outlier## [1] 5.1429durout <- boxplot.stats(pub_mw$dailysleep_hrs_mean)$out
duroutid <- which(pub_mw$dailysleep_hrs_mean %in% c(durout))
pub_mw[duroutid, ]$ELS_ID## [1] 96
## 162 Levels: 1 2 3 5 6 12 13 14 15 16 17 20 21 22 23 24 25 26 29 31 32 33 ... 312pub_mw_rem <-
pub_mw %>%
filter(
!ELS_ID == "96"
)3.1 centering based on this sample
I standardized tanner stage, age at tanner assessment, and bmi within each sex
pub_mw_cent_hrs <-
pub_mw_rem %>%
mutate(
Sex = factor(Sex),
COVID_AHC_MW = factor(COVID_AHC_MW),
mw_dailysleep_sat_z = scale(dailysleep_sat_mean),
mw_dailysleepshrs_z = scale(dailysleep_hrs_mean),
household_income_mw = scale(Household_Income_MW),
sumsev_type_z = scale(sumsev_type),
rr_z = scale(sleep_resp_rate.x)
) %>%
group_by(Sex) %>%
mutate(
tanner_average_for_mw_z = scale(tanner_average_for_mw), # standardize within each sex
tanner_adrenal_mw_z = scale(tanner_adrenal_mw),
tanner_gonadal_mw_z = scale(tanner_gonadal_mw),
bmi_at_mw_z = scale(BMI_MW),
tanner_age_for_mw_z = scale(tanner_age_for_mw)
) %>%
ungroup()3.2 computing relative pubertal stage
rel_pub_mod_hrs <- lm(tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_cent_hrs)
summary(rel_pub_mod_hrs)##
## Call:
## lm(formula = tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_cent_hrs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.400 -0.778 0.171 0.647 1.473
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.000000000000000452 0.087921719844392565 0.00 1
## tanner_age_for_mw_z 0.386878313887527692 0.088732075406578556 4.36 0.00003
##
## (Intercept)
## tanner_age_for_mw_z ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.922 on 108 degrees of freedom
## Multiple R-squared: 0.15, Adjusted R-squared: 0.142
## F-statistic: 19 on 1 and 108 DF, p-value: 0.0000298pub_mw_cent_hrs <-
pub_mw_cent_hrs %>%
add_residuals(rel_pub_mod_hrs, var = "rel_pub_stg_hrs")3.3 associations with hours
testing for potential covariates
cor.test(pub_mw_cent_hrs$mw_dailysleepshrs_z, pub_mw_cent_hrs$bmi_at_mw_z) # not sig##
## Pearson's product-moment correlation
##
## data: pub_mw_cent_hrs$mw_dailysleepshrs_z and pub_mw_cent_hrs$bmi_at_mw_z
## t = -0.788, df = 103, p-value = 0.43
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.26510 0.11601
## sample estimates:
## cor
## -0.077371# r=-0.08045863, p=0.4146
cor.test(pub_mw_cent_hrs$mw_dailysleepshrs_z, pub_mw_cent_hrs$household_income_mw) # not sig##
## Pearson's product-moment correlation
##
## data: pub_mw_cent_hrs$mw_dailysleepshrs_z and pub_mw_cent_hrs$household_income_mw
## t = 1.56, df = 101, p-value = 0.12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.040983 0.337251
## sample estimates:
## cor
## 0.15376# r=0.1545033, p=0.1192
cor.test(pub_mw_cent_hrs$mw_dailysleepshrs_z, pub_mw_cent_hrs$sumsev_type_z) # not sig##
## Pearson's product-moment correlation
##
## data: pub_mw_cent_hrs$mw_dailysleepshrs_z and pub_mw_cent_hrs$sumsev_type_z
## t = -2.77, df = 107, p-value = 0.0067
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.425800 -0.073883
## sample estimates:
## cor
## -0.25839# r=-0.2487005 , p=0.009116
cor.test(pub_mw_cent_hrs$mw_dailysleepshrs_z, pub_mw_cent_hrs$rr_z) # not sig##
## Pearson's product-moment correlation
##
## data: pub_mw_cent_hrs$mw_dailysleepshrs_z and pub_mw_cent_hrs$rr_z
## t = 1.71, df = 108, p-value = 0.09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.025749 0.339209
## sample estimates:
## cor
## 0.16227# r=0.16227 , p=0.09
summary(lm(mw_dailysleepshrs_z ~ Sex, data = pub_mw_cent_hrs)) # not sig##
## Call:
## lm(formula = mw_dailysleepshrs_z ~ Sex, data = pub_mw_cent_hrs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5064 -0.4951 -0.0123 0.7117 1.9186
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.110 0.153 -0.72 0.47
## Sex2 0.181 0.196 0.93 0.36
##
## Residual standard error: 1 on 108 degrees of freedom
## Multiple R-squared: 0.00789, Adjusted R-squared: -0.0013
## F-statistic: 0.858 on 1 and 108 DF, p-value: 0.356# B=0.1915, p=0.329
summary(lm(mw_dailysleepshrs_z ~ COVID_AHC_MW, data = pub_mw_cent_hrs)) # not sig##
## Call:
## lm(formula = mw_dailysleepshrs_z ~ COVID_AHC_MW, data = pub_mw_cent_hrs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4773 -0.5378 0.0513 0.6835 2.0482
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0589 0.1206 -0.49 0.63
## COVID_AHC_MW1 0.1579 0.1975 0.80 0.43
##
## Residual standard error: 1 on 108 degrees of freedom
## Multiple R-squared: 0.00589, Adjusted R-squared: -0.00332
## F-statistic: 0.639 on 1 and 108 DF, p-value: 0.426# B=0.17485, p=0.3784 Tanner Average and Sleep Duration Association
4.1 association between relative pubertal stage score and sleep duration
# relative pubertal stage
pub_mw_cent_hrs_rm <-
pub_mw_cent_hrs %>%
drop_na(sumsev_type)
relpubhrs <- lm(mw_dailysleepshrs_z ~ rel_pub_stg_hrs + sumsev_type_z, data = pub_mw_cent_hrs_rm)
summary(gvlma(rel_pub_mod_hrs))##
## Call:
## lm(formula = tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_cent_hrs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.400 -0.778 0.171 0.647 1.473
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.000000000000000452 0.087921719844392565 0.00 1
## tanner_age_for_mw_z 0.386878313887527692 0.088732075406578556 4.36 0.00003
##
## (Intercept)
## tanner_age_for_mw_z ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.922 on 108 degrees of freedom
## Multiple R-squared: 0.15, Adjusted R-squared: 0.142
## F-statistic: 19 on 1 and 108 DF, p-value: 0.0000298
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = rel_pub_mod_hrs)
##
## Value p-value Decision
## Global Stat 5.695554 0.2231 Assumptions acceptable.
## Skewness 3.960436 0.0466 Assumptions NOT satisfied!
## Kurtosis 1.033864 0.3093 Assumptions acceptable.
## Link Function 0.000902 0.9760 Assumptions acceptable.
## Heteroscedasticity 0.700352 0.4027 Assumptions acceptable.given certain linear relations were violated, testing quadratic term of relative pubertal stage
library(moments)
skewness(pub_mw_cent_hrs_rm$dailysleep_hrs_mean) # -.33929## [1] -0.33929pub_mw_cent_hrs_rm <-
pub_mw_cent_hrs_rm %>%
mutate(
dailysleep_hrs_log = log10(max(dailysleep_hrs_mean + 1) - dailysleep_hrs_mean),
dailysleep_hrs_sqrt = sqrt(max(dailysleep_hrs_mean+1) - dailysleep_hrs_mean)
)
skewness(pub_mw_cent_hrs_rm$dailysleep_hrs_log) # -0.5227 ooo terrible## [1] -0.5227skewness(pub_mw_cent_hrs_rm$dailysleep_hrs_sqrt) # -0.075814 better## [1] -0.075814pub_mw_cent_hrs_rm <- pub_mw_cent_hrs_rm %>% mutate(dailysleep_hrs_sqrt_z = scale(dailysleep_hrs_sqrt))
relpubhrssqrt<- lm(dailysleep_hrs_sqrt_z ~ rel_pub_stg_hrs + sumsev_type_z, data = pub_mw_cent_hrs_rm)
summary(gvlma(relpubhrssqrt))##
## Call:
## lm(formula = dailysleep_hrs_sqrt_z ~ rel_pub_stg_hrs + sumsev_type_z,
## data = pub_mw_cent_hrs_rm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7533 -0.6554 -0.0053 0.6757 2.2534
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0000251 0.0932859 0.00 0.9998
## rel_pub_stg_hrs -0.0075593 0.1019744 -0.07 0.9410
## sumsev_type_z 0.2631559 0.0939661 2.80 0.0061 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.974 on 106 degrees of freedom
## Multiple R-squared: 0.069, Adjusted R-squared: 0.0515
## F-statistic: 3.93 on 2 and 106 DF, p-value: 0.0226
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = relpubhrssqrt)
##
## Value p-value Decision
## Global Stat 2.485152 0.647 Assumptions acceptable.
## Skewness 0.024764 0.875 Assumptions acceptable.
## Kurtosis 0.000371 0.985 Assumptions acceptable.
## Link Function 2.028813 0.154 Assumptions acceptable.
## Heteroscedasticity 0.431205 0.511 Assumptions acceptable.summary(relpubhrssqrt)##
## Call:
## lm(formula = dailysleep_hrs_sqrt_z ~ rel_pub_stg_hrs + sumsev_type_z,
## data = pub_mw_cent_hrs_rm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7533 -0.6554 -0.0053 0.6757 2.2534
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0000251 0.0932859 0.00 0.9998
## rel_pub_stg_hrs -0.0075593 0.1019744 -0.07 0.9410
## sumsev_type_z 0.2631559 0.0939661 2.80 0.0061 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.974 on 106 degrees of freedom
## Multiple R-squared: 0.069, Adjusted R-squared: 0.0515
## F-statistic: 3.93 on 2 and 106 DF, p-value: 0.0226# rel_pub_stg_hrs B=-0.0075593 p=0.9410
# sumsev_type_z B=0.2631559 p=0.0061 **
# Mod R2 = 0.0515, p=0.0226
# CI
ci_relpub = Boot(relpubhrssqrt, coef,
method='case', R=1000)
boot_confint_relpub <- confint(ci_relpub, type='norm')
# rel_pub_stg_hrs -0.22551 0.21573
# sumsev_type_z 0.10225 0.41852
# Bayes Factor
bayes_relpub <- regressionBF(
formula = dailysleep_hrs_sqrt_z ~ rel_pub_stg_hrs + sumsev_type_z, data = pub_mw_cent_hrs_rm
)## Warning: data coerced from tibble to data frame# BF
# [1] rel_pub_stg_hrs : 0.20437 ±0%
# [2] sumsev_type_z : 6.5706 ±0%
hrssumsevint <- summary(relpubhrssqrt)$coefficients[1,1]
sumsvslp <- summary(relpubhrssqrt)$coefficients[3,1]
pub_mw_cent_hrs_rm %>%
ggplot(
aes(
x = sumsev_type_z,
y = dailysleep_hrs_sqrt_z
)
) +
geom_point(size = 2, alpha = .5) +
geom_abline(
intercept = hrssumsevint,
slope = sumsvslp,
size=2
) +
theme_classic() +
labs(
x = "Severity of ELS (z-scored)",
y = "Sleep Duration (subjective)",
title = "Association between severity of ELS and Sleep Duration"
) +
theme(
axis.text = element_text(size = 10, angle = 30, hjust = 1),
axis.title = element_text(size = 10),
plot.title = element_text(size = 10)
)
ggsave("els_sleep_dur_subj.png")## Saving 7 x 5 in image5 Circadian preference
pub_mw_cmep <-
pub_mw %>%
drop_na(cmep_total_mw)
summary(pub_mw_cmep$cmep_total_mw)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 13.0 21.0 26.0 25.4 29.0 39.0 # Min. 1st Qu. Median Mean 3rd Qu. Max.
# 13.00 21.00 26.00 25.43 29.00 39.00
sd(pub_mw_cmep$cmep_total_mw) # 5.558876## [1] 5.55895.1 outliers - circadian preference
boxplot.stats(pub_mw_cmep$cmep_total_mw)$out # 0 outliers## numeric(0)cor.test(pub_mw_cmep$cmep_total_mw, pub_mw_cmep$dailysleep_sat_mean) # no association between sleep satisfaction and cmep##
## Pearson's product-moment correlation
##
## data: pub_mw_cmep$cmep_total_mw and pub_mw_cmep$dailysleep_sat_mean
## t = 3.98, df = 103, p-value = 0.00013
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.18605 0.52002
## sample estimates:
## cor
## 0.36471# r=0.36471, p=0.00013
cor.test(pub_mw_cmep$cmep_total_mw, pub_mw_cmep$dailysleep_hrs_mean)##
## Pearson's product-moment correlation
##
## data: pub_mw_cmep$cmep_total_mw and pub_mw_cmep$dailysleep_hrs_mean
## t = 2.85, df = 103, p-value = 0.0052
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.083404 0.439595
## sample estimates:
## cor
## 0.27074# r=0.27074, p=0.00525.2 centering based on this sample
I standardized tanner stage, age at tanner assessment, and bmi within each sex
pub_mw_cmep_cent <-
pub_mw_cmep %>%
mutate(
Sex = factor(Sex),
COVID_AHC_MW = factor(COVID_AHC_MW),
mw_dailysleep_sat_z = scale(dailysleep_sat_mean),
mw_dailysleepshrs_z = scale(dailysleep_hrs_mean),
household_income_mw = scale(Household_Income_MW),
cmep_total_mw_z = scale(cmep_total_mw),
sumsev_type_z = scale(sumsev_type),
rr_z = scale(sleep_resp_rate.x)
) %>%
group_by(Sex) %>%
mutate(
tanner_average_for_mw_z = scale(tanner_average_for_mw), # standardize within each sex
tanner_adrenal_mw_z = scale(tanner_adrenal_mw),
tanner_gonadal_mw_z = scale(tanner_gonadal_mw),
bmi_at_mw_z = scale(BMI_MW),
tanner_age_for_mw_z = scale(tanner_age_for_mw)
) %>%
ungroup()5.3 computing relative pubertal stage
rel_pub_modcmep <- lm(tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_cmep_cent)
summary(rel_pub_modcmep)##
## Call:
## lm(formula = tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_cmep_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.330 -0.825 0.171 0.656 1.436
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.000000000000000358 0.090252870398000395 0.00
## tanner_age_for_mw_z 0.380412980084517016 0.091124899061470935 4.17
## Pr(>|t|)
## (Intercept) 1
## tanner_age_for_mw_z 0.000063 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.925 on 103 degrees of freedom
## Multiple R-squared: 0.145, Adjusted R-squared: 0.136
## F-statistic: 17.4 on 1 and 103 DF, p-value: 0.0000625pub_mw_cmep_cent <-
pub_mw_cmep_cent %>%
add_residuals(rel_pub_modcmep, var = "rel_pub_stg_cmep")5.4 associations with possible covariates
cor.test(pub_mw_cmep_cent$cmep_total_mw_z, pub_mw_cmep_cent$bmi_at_mw_z) # not sig##
## Pearson's product-moment correlation
##
## data: pub_mw_cmep_cent$cmep_total_mw_z and pub_mw_cmep_cent$bmi_at_mw_z
## t = 0.691, df = 98, p-value = 0.49
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.12852 0.26248
## sample estimates:
## cor
## 0.069653# r=0.069653, p=0.49
cor.test(pub_mw_cmep_cent$cmep_total_mw_z, pub_mw_cmep_cent$household_income_mw) # not sig##
## Pearson's product-moment correlation
##
## data: pub_mw_cmep_cent$cmep_total_mw_z and pub_mw_cmep_cent$household_income_mw
## t = 0.195, df = 97, p-value = 0.85
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.17829 0.21639
## sample estimates:
## cor
## 0.01982# r=0.01982, p=0.85
cor.test(pub_mw_cmep_cent$cmep_total_mw_z, pub_mw_cmep_cent$sumsev_type_z) # not sig##
## Pearson's product-moment correlation
##
## data: pub_mw_cmep_cent$cmep_total_mw_z and pub_mw_cmep_cent$sumsev_type_z
## t = -0.0915, df = 103, p-value = 0.93
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.20034 0.18296
## sample estimates:
## cor
## -0.0090194# r=-0.0090194, p=0.93
cor.test(pub_mw_cmep_cent$cmep_total_mw_z, pub_mw_cmep_cent$rr_z) # not sig##
## Pearson's product-moment correlation
##
## data: pub_mw_cmep_cent$cmep_total_mw_z and pub_mw_cmep_cent$rr_z
## t = 2.12, df = 103, p-value = 0.037
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.013102 0.381004
## sample estimates:
## cor
## 0.20425# r=0.20425, p=0.037
summary(lm(cmep_total_mw_z ~ Sex, data = pub_mw_cmep_cent)) # moderate association ##
## Call:
## lm(formula = cmep_total_mw_z ~ Sex, data = pub_mw_cmep_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1383 -0.6991 0.0204 0.7400 2.1792
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.262 0.148 1.78 0.079 .
## Sex2 -0.451 0.194 -2.33 0.022 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.979 on 103 degrees of freedom
## Multiple R-squared: 0.0501, Adjusted R-squared: 0.0409
## F-statistic: 5.43 on 1 and 103 DF, p-value: 0.0217# B=-0.451, p=0.022
summary(lm(cmep_total_mw_z ~ COVID_AHC_MW, data = pub_mw_cmep_cent)) # not sig##
## Call:
## lm(formula = cmep_total_mw_z ~ COVID_AHC_MW, data = pub_mw_cmep_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0757 -0.6365 0.0082 0.7277 2.3468
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0946 0.1227 0.77 0.44
## COVID_AHC_MW1 -0.2547 0.2014 -1.26 0.21
##
## Residual standard error: 0.997 on 103 degrees of freedom
## Multiple R-squared: 0.0153, Adjusted R-squared: 0.00574
## F-statistic: 1.6 on 1 and 103 DF, p-value: 0.209# B=-0.2547, p=0.21cmepmod <- lm(cmep_total_mw_z ~ rel_pub_stg_cmep + Sex + rr_z + mw_dailysleepshrs_z + mw_dailysleep_sat_z, data = pub_mw_cmep_cent)
summary(gvlma(cmepmod)) # all checks out##
## Call:
## lm(formula = cmep_total_mw_z ~ rel_pub_stg_cmep + Sex + rr_z +
## mw_dailysleepshrs_z + mw_dailysleep_sat_z, data = pub_mw_cmep_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1074 -0.6932 -0.0427 0.5850 2.1266
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2774 0.1364 2.03 0.0447 *
## rel_pub_stg_cmep 0.0879 0.1021 0.86 0.3916
## Sex2 -0.4775 0.1806 -2.64 0.0095 **
## rr_z 0.2232 0.0920 2.43 0.0171 *
## mw_dailysleepshrs_z 0.1027 0.1059 0.97 0.3343
## mw_dailysleep_sat_z 0.3192 0.1071 2.98 0.0036 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.894 on 99 degrees of freedom
## Multiple R-squared: 0.24, Adjusted R-squared: 0.202
## F-statistic: 6.25 on 5 and 99 DF, p-value: 0.0000441
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = cmepmod)
##
## Value p-value Decision
## Global Stat 3.739 0.443 Assumptions acceptable.
## Skewness 0.167 0.683 Assumptions acceptable.
## Kurtosis 1.189 0.275 Assumptions acceptable.
## Link Function 0.309 0.579 Assumptions acceptable.
## Heteroscedasticity 2.074 0.150 Assumptions acceptable.car::vif(cmepmod) # all below 5## rel_pub_stg_cmep Sex rr_z mw_dailysleepshrs_z
## 1.1511 1.0440 1.1026 1.4598
## mw_dailysleep_sat_z
## 1.4938summary(cmepmod)##
## Call:
## lm(formula = cmep_total_mw_z ~ rel_pub_stg_cmep + Sex + rr_z +
## mw_dailysleepshrs_z + mw_dailysleep_sat_z, data = pub_mw_cmep_cent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1074 -0.6932 -0.0427 0.5850 2.1266
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2774 0.1364 2.03 0.0447 *
## rel_pub_stg_cmep 0.0879 0.1021 0.86 0.3916
## Sex2 -0.4775 0.1806 -2.64 0.0095 **
## rr_z 0.2232 0.0920 2.43 0.0171 *
## mw_dailysleepshrs_z 0.1027 0.1059 0.97 0.3343
## mw_dailysleep_sat_z 0.3192 0.1071 2.98 0.0036 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.894 on 99 degrees of freedom
## Multiple R-squared: 0.24, Adjusted R-squared: 0.202
## F-statistic: 6.25 on 5 and 99 DF, p-value: 0.0000441# rel_pub_stg_cmep B=0.0879, p=0.3916
# Sex2 B=-0.4775, p=0.0095
# rr_z B = 0.2232, p=0.0171
# mw_dailysleepshrs_z B=0.1027, p=0.3343
# mw_dailysleep_sat_z B=0.3192, p=0.0036
# overall R2=0.202, p<0.001
cmepboot = Boot(cmepmod, coef,
method='case', R=1000)
# bootstrap standard confidence interval
cmepbootconfint <- confint(cmepboot, type='norm')
# rel_pub_stg_cmep -0.115325 0.30192
# Sex2 -0.790891 -0.13460
# rr_z 0.037564 0.43326
# mw_dailysleepshrs_z -0.080480 0.27529
# mw_dailysleep_sat_z 0.140412 0.50786
bayes <- generalTestBF(
formula = cmep_total_mw_z ~ rel_pub_stg_cmep + Sex + rr_z + mw_dailysleepshrs_z + mw_dailysleep_sat_z, data = pub_mw_cmep_cent
)## Warning: data coerced from tibble to data frame# --------------
# [1] rel_pub_stg_cmep: 0.21759
# [2] Sex: 2.2676
# [3] rr_z: 1.4939
# [4] mw_dailysleepshrs_z: 7.2318
# [8] mw_dailysleep_sat_z: 173.346 correlations among sleep measures
6.1 correlations
library(modelr)
library(corrr)
pub_mw_cmep_cent %>%
dplyr::select(
cmep_total_mw,
dailysleep_sat_mean,
dailysleep_hrs_mean
) %>%
correlate() %>%
fashion()##
## Correlation method: 'pearson'
## Missing treated using: 'pairwise.complete.obs'## term cmep_total_mw dailysleep_sat_mean dailysleep_hrs_mean
## 1 cmep_total_mw .36 .27
## 2 dailysleep_sat_mean .36 .48
## 3 dailysleep_hrs_mean .27 .48cor.test(pub_mw_cmep_cent$cmep_total_mw, pub_mw_cmep_cent$dailysleep_sat_mean) # r = 0.3619217, p<0.001##
## Pearson's product-moment correlation
##
## data: pub_mw_cmep_cent$cmep_total_mw and pub_mw_cmep_cent$dailysleep_sat_mean
## t = 3.98, df = 103, p-value = 0.00013
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.18605 0.52002
## sample estimates:
## cor
## 0.36471cor.test(pub_mw_cmep_cent$cmep_total_mw, pub_mw_cmep_cent$dailysleep_hrs_mean) # r = 0.27074, p=0.0052##
## Pearson's product-moment correlation
##
## data: pub_mw_cmep_cent$cmep_total_mw and pub_mw_cmep_cent$dailysleep_hrs_mean
## t = 2.85, df = 103, p-value = 0.0052
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.083404 0.439595
## sample estimates:
## cor
## 0.27074corr_mw_sleep <-
pub_mw_cmep_cent %>%
dplyr::select(
cmep_total_mw,
dailysleep_sat_mean,
dailysleep_hrs_mean
) %>%
rename(
`Circadian Preference` = cmep_total_mw,
`Sleep Satisfaction` = dailysleep_sat_mean,
`Sleep Duration` = dailysleep_hrs_mean
) %>%
corrr::correlate() %>%
corrr::shave(upper = FALSE)##
## Correlation method: 'pearson'
## Missing treated using: 'pairwise.complete.obs'corr_mw_sleep## # A tibble: 3 x 4
## term `Circadian Preferenc… `Sleep Satisfactio… `Sleep Duration`
## <chr> <dbl> <dbl> <dbl>
## 1 Circadian Preferen… NA 0.365 0.271
## 2 Sleep Satisfaction NA NA 0.478
## 3 Sleep Duration NA NA NAcorr_mw_sleep_long <-
corr_mw_sleep %>%
pivot_longer(
cols = -term,
names_to = "colname",
values_to = "corr"
) %>%
mutate(
rowname = fct_inorder(term),
colname = fct_inorder(colname)
)
corr_mw_sleep_long## # A tibble: 9 x 4
## term colname corr rowname
## <chr> <fct> <dbl> <fct>
## 1 Circadian Preference Circadian Preference NA Circadian Preference
## 2 Circadian Preference Sleep Satisfaction 0.365 Circadian Preference
## 3 Circadian Preference Sleep Duration 0.271 Circadian Preference
## 4 Sleep Satisfaction Circadian Preference NA Sleep Satisfaction
## 5 Sleep Satisfaction Sleep Satisfaction NA Sleep Satisfaction
## 6 Sleep Satisfaction Sleep Duration 0.478 Sleep Satisfaction
## 7 Sleep Duration Circadian Preference NA Sleep Duration
## 8 Sleep Duration Sleep Satisfaction NA Sleep Duration
## 9 Sleep Duration Sleep Duration NA Sleep Duration# plotting correlations
corrplotsleep <-
corr_mw_sleep_long %>%
ggplot(
aes(
x = rowname,
y = fct_rev(colname),
fill = corr
)
) +
geom_tile() +
geom_text(
aes(
label = round(corr, 2)
)
) +
coord_fixed(expand = FALSE) +
scale_fill_distiller(
palette = "RdBu",
na.value = "white",
direction = 1,
limits = c(-1,1)
) +
labs(
x = "",
y = ""
) +
theme(
axis.text = element_text(size = 10, angle = 30, hjust = 1)
)
corrplotsleep## Warning: Removed 6 rows containing missing values (geom_text).
ggsave("corrplot_sleep_mw.png", corrplotsleep)## Saving 7 x 5 in image## Warning: Removed 6 rows containing missing values (geom_text).emasleeplongfp <- "~/Box/Mooddata_Coordinating/1_Lab_Coordinating/Users/JackieSchwartz/Dissertation/0_MW_Act_Demo_Descriptives/MW_daily_sleep_data_longform.csv"
emasleeplong <-
read_csv(emasleeplongfp) %>%
mutate(ELS_ID = factor(ELS_ID))
emasleeplong_2wkpub <-
left_join(
emasleeplong,
pub_mw_cent,
by = "ELS_ID"
) %>%
drop_na(tanner_average_for_mw)
emasleeplong_2wkpub %>%
distinct(ELS_ID)## # A tibble: 111 x 1
## ELS_ID
## <fct>
## 1 2
## 2 69
## 3 20
## 4 14
## 5 75
## 6 50
## 7 58
## 8 59
## 9 95
## 10 77
## # … with 101 more rowswrite_csv(emasleeplong_2wkpub, "emasleeplong_2wk_final.csv")6.1.1 distribution of weekday vs. weekend and association with sleep sat and hrs
emasleeplong_2wkpub %>%
ggplot(
aes(x = week, y = DailySleep_satisfaction, fill = week)
) +
geom_violin(alpha=0.5, color= "black") +
geom_boxplot(width=0.1, color = "grey", alpha=0.5) +
scale_fill_manual(values = c("#FFA07A","#6B8E23")) +
theme_classic() +
labs(x = "Weekday or Weekend", y = "Sleep Satisfaction")
ggsave("wkdayvsweeknd_sleepsat.png", width = 7, height = 6)
summary(lmer(scale(DailySleep_satisfaction) ~ scale(dayorder) +factor(week) + (1 + scale(dayorder)|ELS_ID), data = emasleeplong_2wkpub))## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: scale(DailySleep_satisfaction) ~ scale(dayorder) + factor(week) +
## (1 + scale(dayorder) | ELS_ID)
## Data: emasleeplong_2wkpub
##
## REML criterion at convergence: 2485.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.209 -0.558 0.073 0.577 3.172
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ELS_ID (Intercept) 0.4183 0.647
## scale(dayorder) 0.0243 0.156 0.18
## Residual 0.5656 0.752
## Number of obs: 985, groups: ELS_ID, 111
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.05185 0.06884 117.24854 -0.75 0.45285
## scale(dayorder) -0.00816 0.03271 55.27765 -0.25 0.80405
## factor(week)weekend 0.19116 0.05786 873.30592 3.30 0.00099 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(d)
## scal(dyrdr) 0.163
## fctr(wk)wkn -0.192 0.058emasleeplong_2wkpub %>%
ggplot(
aes(x = week, y = as.numeric(DailySleep_hrs_rec), fill = week)
) +
geom_violin(width=.8, alpha=0.5, color= "black") +
geom_boxplot(width=0.03, color = "grey", alpha=0.5) +
scale_fill_manual(values = c("#FFA07A","#6B8E23")) +
theme_classic() +
labs(x = "Weekday or Weekend", y = "Subjective Sleep Duration")
ggsave("wkdayvsweeknd_sleephrs.png", width = 7, height = 6)
summary(lmer(scale(DailySleep_hrs_rec) ~ scale(dayorder) + factor(week) + (1 + scale(dayorder)|ELS_ID), data = emasleeplong_2wkpub))## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: scale(DailySleep_hrs_rec) ~ scale(dayorder) + factor(week) +
## (1 + scale(dayorder) | ELS_ID)
## Data: emasleeplong_2wkpub
##
## REML criterion at convergence: 2508.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9882 -0.5592 0.0381 0.6480 2.6954
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ELS_ID (Intercept) 0.400 0.632
## scale(dayorder) 0.010 0.100 -0.29
## Residual 0.591 0.769
## Number of obs: 985, groups: ELS_ID, 111
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.0695 0.0674 116.6825 -1.03 0.304
## scale(dayorder) 0.0346 0.0292 57.0247 1.19 0.241
## factor(week)weekend 0.1495 0.0588 885.0031 2.54 0.011 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(d)
## scal(dyrdr) -0.024
## fctr(wk)wkn -0.203 0.055sleep_week <-
emasleeplong_2wkpub %>%
group_by(ELS_ID, week) %>%
summarise(
n = n()
)## `summarise()` has grouped output by 'ELS_ID'. You can override using the `.groups` argument.emasleep_filter_spread <-
emasleeplong_2wkpub %>%
spread(week, DailySleep_satisfaction)
emasleeplong_2wkpub %>%
ggplot(
aes(
x = dayorder,
y = DailySleep_satisfaction,
color = week
)
) +
geom_line(aes(group = ELS_ID, color = week), alpha = .3) +
geom_smooth(method = "loess", se=FALSE, size = 2) +
theme_classic() +
scale_x_continuous(
name = "Day",
limits = c(1, 14),
breaks = seq(1, 14, 1)
) +
scale_y_continuous(
name = "Sleep Satisfaction",
limits = c(1, 100),
breaks = seq(1, 100, 10)
) +
scale_color_manual(values = c("#FFA07A","#6B8E23"))## `geom_smooth()` using formula 'y ~ x'
ggsave("sleepsat_byDay_byWkdayWkend.png", width = 7, height = 6)## `geom_smooth()` using formula 'y ~ x'emasleeplong_2wkpub %>%
ggplot(
aes(
x = dayorder,
y = DailySleep_hrs_rec,
color = week
)
) +
geom_line(aes(group = ELS_ID, color = week), alpha = .3) +
geom_smooth(method = "loess", se=FALSE, size = 2) +
theme_classic() +
scale_x_continuous(
name = "Day",
limits = c(1, 14),
breaks = seq(1, 14, 1)
) +
scale_y_continuous(
name = "Sleep Satisfaction",
limits = c(5, 9),
breaks = seq(5, 9, 1)
) +
scale_color_manual(values = c("#FFA07A","#6B8E23"))## `geom_smooth()` using formula 'y ~ x'
ggsave("sleephrs_byDay_byWkdayWkend.png", width = 7, height = 6)## `geom_smooth()` using formula 'y ~ x'